function beta = BsplineCDFLsqFit(samp,knotvec,EvalPts,SampWts)

if sum(samp<knotvec(1))>0 || sum(samp>knotvec(end))>0
    error('BsplineCDFLsqFit:BadKnotVec','THE KNOT VECTOR MUST SPAN THE SAMPLE SPACE.')
end

if length(EvalPts) < length(knotvec)
    error('BsplineCDFLsqFit:NotIdentified','TOO FEW ECDF EvalPts: CANNOT ESTIMATE A NON-IDENTIFIED MODEL.')
end

if nargin<4
    SampWts = ones(length(samp),1)/length(samp) ;
else
    SampWts = SampWts/sum(SampWts) ;
end

Fhat = zeros(size(EvalPts));
for i=1:length(EvalPts)
    Fhat(i) = sum( SampWts.*(samp<=EvalPts(i)) ) ;
end

[C] = BsplineBasis3(knotvec,EvalPts,'off');
A = [];
b = [];
K = length(knotvec)-1;
for i=1:K+2  %%%%This loop will give us the monotonicity constraints.
    A = [A; zeros(1,i-1),1,-1,zeros(1,(K+3)-2-(i-1))];
    b = [b; 0];
end

Aeq = BsplineBasis3(knotvec,[knotvec(1);knotvec(end)],'off');
beq = [0,1];
lb = zeros(K+3,1);
ub = ones(K+3,1);

options = optimoptions('lsqlin','Algorithm','interior-point','Display','none');

beta = lsqlin(C,Fhat,A,b,Aeq,beq,lb,ub,[],options);


